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Abstract 

A DNS  code  based  on  a spectral-element  method  is  developed.  The  code 
is  applied  to  calculate  Eulerian  space-time  velocity  correlation  functions 
and  to  study  the  possibility  to  describe  the  dispersion  of  particles  by  a 
Langevin  equation  for  particle  velocity.  The  preliminary  results  indicate 
that  it  is  possible  to  derive  the  damping  coefficients  from  the  correlation 
functions. 

1.  Introduction 

Particle  dispersion  can  be  described  by  a Langevin  equation,  but  there  is 
no  unique  way  to  determine  the  coefficients  in  this  equation  for  inhomoge- 
neous flow.  A recently  developed  theory  (Brouwers  [1])  predicts  that  these 
coefficients  are  related  to  the  Eulerian  space-time  correlation  function  in 
a moving  frame.  The  frame  is  moved  in  time  with  the  local  average  fluid 
velocity. 

The  outline  of  this  paper  is  as  follows.  In  section  2 the  numerical  methods 
for  the  continuum  phase  and  particles  are  considered.  Next,  in  section  3 
the  stochastic  theory  is  explained.  Then  some  preliminary  results  will  be 
presented  and  discussed  in  section  4.  Finally,  the  conclusions  are  made. 

2.  Numerical  Method 

2.1.  CONTINUUM  PHASE 

In  this  section  the  numerical  method  used  for  DNS  of  turbulent  flow  in  a 
cylindrical  pipe  will  be  described.  The  equations  to  be  solved  are  the  Navier- 
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Stokes  equations  and  continuity  equation  for  incompressible  flow.  Because 
of  the  cylindrical  geometry  the  choice  for  cylindrical  coordinates  and  cylin- 
drical velocity  components  is  natural.  Used  are  r,  cp  and  z for  the  radial, 
tangential  and  axial  coordinates  and  ur , and  uz  for  the  corresponding 
velocity  components.  Then  the  governing  equations  can  be  written  in  the 
form: 


r or  r o<p  oz 


(1) 


< 


dur  dP 

— + u )<i>uz  - uzu^  + = 


duj, 

dt 


-f-  LOzUj-  — 0JTUz  -f-  — 


1 dP 


r d<j> 


2 ( \ 2 ^ 
Re\AUr~72'P2^) 


_2_ 

Re 


(Au<t>  - ^ + 


2 dur\ 
r2  d<p  ) 


diiz  dP  2 

[ -gj-  + UrU^  - WtfiUr  + ~fa  = ~ReAUz  + * 


(2) 


Here  t denotes  time,  ajr,  and  u)z  are  the  cylindrical  components  of  the 
vorticity,  A is  the  Laplace  operator  and  the  total  pressure  P is  defined  as 
P = p + u2/2,  where  p is  the  static  pressure.  The  forcing  term  / is  adjusted 
each  time  step  to  ensure  a constant  mass  flow  through  the  pipe  demanding 


/ = - 


_4_  dv^ 
Re  dr 


r=R 


The  equations  have  been  non-dimensionalized  by  the  radius  of  the  pipe  R, 
the  kinematic  viscosity  u and  the  bulk  velocity  u b ■ The  Reynolds  number 
is  defined  as  Re  = uqD/v  where  D is  the  diameter  of  the  pipe. 

A finite  part  of  a cylindrical  pipe  of  length  L is  considered,  using  periodic 
boundary  conditions  in  the  axial  direction.  Combined  with  the  natural  pe- 
riodicity in  the  tangential  direction  there  are  two  periodic  directions  and 
the  choice  for  a spectral  method  is  obvious.  In  the  radial  direction  an  ex- 
pansion based  on  Chebyshev  polynomials  is  used.  However,  the  distribution 
of  the  Gauss-Lobatto  collocation  points  leads  to  very  small  cells  near  the 
pipe  axis  and  thus  necessitates  the  use  of  very  small  time  steps.  Therefore, 
the  radial  direction  is  divided  into  several  elements  and  in  each  element 
an  expansion  into  Chebyshev  polynomials  is  adopted.  At  the  interfaces  be- 
tween the  elements  the  solution  is  required  to  be  Cl . Each  component  of 
the  solution  is  thus  expanded  as: 


M^/2—l  Mz/  2-1 

ti(r,  (p,  z)  = E E u^,kz(r)  exp(ik<p4>  + 2irikzz/L). 

2+1  fc2=— AU/2+1 
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In  this  way  a hybrid  method  appears:  Fourier-Galerkin  in  the  two  periodic 
directions  and  Chebyshev-collocation  in  the  radial  direction.  Derivatives  in 
the  periodic  directions  can  easily  be  calculated  in  spectral  space,  whereas 
derivatives  with  respect  to  r follow  from  the  Chebyshev  derivative  matrix, 
see  [5].  The  reduction  on  the  time  step  is  alleviated  further  by  applying 
’mode-reduction’  in  the  central  element,  i.e.  the  element  near  the  pipe-axis. 
Here,  only  half  of  the  available  modes  are  solved,  see  e.g.  [6].  For  these 
modes  the  boundary  conditions  for  the  velocity  and  pressure  at  r = 0 are 
translated  to  the  element  interface. 

The  Navier-Stokes  equations  are  integrated  in  time  using  a time-splitting 
method  by  [4].  We  write  the  Navier-Stokes  equations  schematically  as 

— + N(u)  + VP  = L(u)  + / , 

where  N denotes  the  nonlinear  terms  on  the  left-hand  side  of  (2),  L the 
viscous  terms  on  the  right-hand  side  and  / the  forcing  term.  A second-order 
accurate  time-splitting  method  with  constant  time  step  At  is  then  given  by 

' un+ 1/3  = 2 un  - in”-1  - At(2N  (un)  - N{un-'))  + Atf 

a pn+i  _ J_y  . ,,n+l/3 

At  v u (‘V \ 

u n+2/3  _ un+ 1/3  _ A/VPn+1  { 

h un+x  = §un+2/3  + §A tL{un+x). 

In  these  formulas  the  superscript  denotes  the  time  level.  In  the  first  step  the 
nonlinear  terms  are  treated  explicitly.  The  products  of  vorticity  and  velocity 
are  calculated  in  a pseudo-spectral  way,  where  fast-fourier  transforms  are 
used  to  transform  from  spectral  to  physical  space.  Note  that  the  nonlinear 
terms  have  to  be  calculated  only  once  per  time  step.  The  latest  result  is 
stored  and  used  again  in  the  next  time  step.  Aliasing  is  prevented  by  the 
3/2-rule. 

In  the  second  step  a Poisson  equation  is  solved  for  the  pressure  to  ensure 
that  the  solution  after  the  third  step  satisfies  the  continuity  equation.  A 
point  that  requires  special  attention  in  this  step  is  the  boundary  condition 
for  the  pressure  at  the  wall  of  the  pipe.  Following  [4]  we  employ  at  r = R 

dPn+1  2 /l^y  d2v%+l\ 

dr  Re  yr  drdcj)  drdz  J ’ 

where  the  solution  at  the  new  time  level  is  found  by  linear  extrapolation 
from  the  solution  at  the  two  previous  time  levels.  In  order  to  obtain  a unique 
solution  this  boundary  condition  cannot  be  used  for  the  {k^^kz)  — (0,0) 
mode.  Instead  the  mean  pressure  at  the  wall  of  the  pipe  is  prescribed. 
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In  the  final  step  the  linear,  viscous  terms  in  the  Navier-Stokes  equations 
are  treated  implicitly.  In  cylindrical  coordinates  the  viscous  parts  of  the 
equations  for  ur  and  are  coupled,  but  they  can  be  decoupled  by  intro- 
ducing u±  = ur  ± iu<p.  Moreover,  the  equations  for  the  Fourier  modes  are 
completely  decoupled,  so  that  a one-dimensional  equation  for  each  Fourier 
mode  and  each  velocity  component  results.  These  one-dimensional  problems 
are  solved  by  a direct  method.  At  the  wall  of  the  pipe  a no-slip  boundary 
condition  is  applied.  At  the  pipe  axis  the  boundary  conditions  follow  from 
the  requirement  that  the  Cartesian  velocity  components  are  regular.  This 
results  in  uz  — 0 for  k,p  ^ 0,  duz/dr  = 0 for  k^  = 0,  u±  = 0 for  k $ / Tl 
and  du±/dr  = 0 for  k $ = =Fl. 

The  calculation  is  started  from  a random  field  superposed  on  an  approxi- 
mate mean  field  with  the  axial  velocity  component  given  by  a logarithmic 
velocity  profile.  The  random  field  is  chosen  in  such  a way  that  it  satisfies  the 
continuity  equation  and  that  the  lowest  four  Fourier  modes  in  both  periodic 
directions  are  unequal  to  zero.  In  the  first  time  step  scheme  (3)  cannot  be 
applied  since  only  one  field  is  available.  A first-order  time-splitting  scheme 
is  used  instead.  After  a large  number  of  time  steps  a state  of  fully-developed 
turbulence  is  reached.  From  that  time  onwards  averaged  flow  quantities  can 
be  calculated. 

2.2.  FLUID  PARTICLES 

Fluid  particles,  e.g.  particles  without  inertia,  are  tracked  using  this  code. 
The  differential  equation  for  the  position  of  such  a particle  is  given  by 

dxi  . . / . , 

— = u{xht),  (4) 

where  Xi  is  the  particle  position  of  particle  i and  u are  the  three  components 
of  the  particle  velocity  derived  at  the  particle  position.  There  are  several 
methods  to  obtain  the  particle  velocity  in  the  interior  of  a grid  cell.  A very 
accurate  method  is  direct  summation  of  the  Fourier  series  expansion,  but 
this  is  prohibitively  expensive.  Linear  interpolation  gives  reasonable  results 
for  the  velocity,  but  derivatives  of  the  velocity  will  be  very  inaccurate. 
Alternatively,  Hermite  interpolation  in  the  radial  direction  and  a fourth 
order  accurate  Lagrange  interpolation  scheme  in  the  periodic  directions  are 
used,  see  [7].  This  is  done  in  physical  space,  since  the  solution  is  known  there 
each  time  step  to  calculate  the  nonlinear  term.  An  Euler-forward  or  2-stage 
Runge-Kutta  method  is  used  to  solve  eq.  (4)  numerically  and  obtain  the 
new  particle  position.  If  a particle  leaves  the  computational  domain  through 
the  axial  boundary  the  periodicity  of  the  velocity  is  used  to  determine  its 
velocity.  If  a particle  leaves  the  pipe  at  the  outer  wall  it  is  reintroduced, 
using  an  elastic  collision  with  the  wall. 
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3.  Stochastic  Approach 

3.1.  DNS 

With  this  code  Eulerian  space-time  velocity  correlation  functions  are  cal- 
culated. These  are  defined  as 

Rii(r , r)  = «{(r,  <f>,  z,  *)u'(r,  t/>,  z + «j(r)r,  t + r),  (5) 

where  the  overhead  bar  means  time-  or  ensemble-average,  «,•  are  the  fluctu- 
ating parts  of  the  three  velocity  components  and  r is  the  time  separation. 
Since  4>  and  z are  coordinates  in  homogeneous  directions  and  the  turbu- 
lence is  stationary,  the  correlation  function  depends  only  on  r and  r.  As 
can  be  seen  from  eq.  (5)  the  frame  is  moved  in  streamwise  direction  with 
the  local  average  axial  velocity.  The  Fourier  transform  of  eq.  (5)  gives  the 
spectral  density  function  for  the  velocity. 

3.2.  LANGEVIN  MODEL 

The  Langevin  equation  for  particle  velocity  can  be  written  as  (see  e.g. 
Borgas  [2]  and  Wilson  &:  Sawford  [3]) 

v'  + r~  V = {C0e)1/2wi{t),  (6) 

where  Co  is  the  universal  Kolomorov  constant,  e the  mean  rate  of  energy 
dissipation  ,w\  ( t ) a white- noise  term  and  rs  the  correlation  time.  The  damp- 
ing coefficient  in  the  Langevin  equation  corresponds  to  the  correlation  time 
ts  of  the  Eulerian  correlation  function.  The  determination  of  ts  from  the 
correlation  function  is  not  a trivial  process,  since  several  time  scales  deter- 
mine the  shape  of  this  function.  The  smallest  scales  cause  the  zero-slope  at 
r = 0.  For  high  Re-numbers  this  effect  decreases  and  the  correlation  func- 
tion could  be  described  by  a single  exponential  function  from  which  the 
correlation  time  can  easily  be  determined.  However,  in  our  simulations  of 
moderate  Re-numbers  a separation  of  time  scales  is  used.  This  enables  us 
to  quantify  the  small  scale  effects. 

To  model  the  second  (small)  time  scale,  smoothing  of  white-noise  is  applied 
(see  e.g.  Stratonovic  [8]),  i.e.  the  rhs  of  eq.  (6)  is  replaced  by  the  random 
process  ((t)  which  obeys  the  equation 

rni  + ^(C0e)^2w2(t),  (7) 

where  tv  is  the  Kolmogorov  time-scale  and  w2{t)  another  white- noise  term. 
The  resulting  spectral  density  of  the  particle  velocity  is 

C0e 

(1  + T2u>2)(rf2  + u2) 


SvH  = 


(8) 
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The  white-noise  from  the  original  Langevin  equation  (eq.  6)  is  no  longer 
lumped  into  a delta-correlated  process.  By  decreasing  the  order  of  the 
white- noise,  viz.  applying  eq.  (7),  the  real  turbulent  process  is  more  ac- 
curately described.  The  process  is  divided  into  2 time  scales.  If  Re  -4  oo, 
7V,  -4  0 and  the  original  Langevin  equation  is  obtained  again.  In  real  flows 
the  Reynolds  number  is  never  infinite,  so  smoothing  the  white-noise  is 
a more  realistic  approach.  The  corresponding  correlation  function  is  the 
Fourier  transform  of  eq.  (8)  and  reads 


t;(0)w(<)  = -—IE  (e  T*  ~^e  Tr> ) > 


(9) 


where  the  standard  deviation  av  is  given  by 


2 _ lr  er8 
V 2C°(1  + ^) 


(10) 


The  spectral  density  function  for  the  acceleration,  Sa  (w),  can  be  obtained 
by  multiplying  eq.  (8)  with  u2. 


4.  Results 

The  investigated  Reynolds  number  equals  5300  based  on  the  bulk  velocity. 
The  number  of  Fourier-modes  in  tangential  and  axial  directions  equal  resp. 
128  and  128.  5 Chebyshev  elements  are  used  which  contain  28,  28,  21,  28 
and  7 collocation  points  resp.,  where  the  element  near  the  pipe  axis  has 
7 collocation  points.  The  code  has  been  extensively  tested  and  compared 
with  DNS  and  experimental  results  of  others  ([6],  [9],  [10]  and  [11]).  The 
investigated  results,  e.g.  rms  of  velocity  fluctuations,  pressure-  and  vorticity 
fluctuations,  2-point  velocity  spectra  and  the  components  of  the  turbulent 
kinetic  energy  equation,  are  in  good  agreement.  These  quantities  are  merely 
used  to  check  our  code  and  hence  the  results  are  not  displayed  here. 

Some  preliminary  results  on  particle  dispersion  will  be  discussed  now.  Spe- 
cial emphasis  is  put  on  the  tangential  velocity  component.  Due  to  the  sym- 
metry of  the  flow,  the  Langevin  for  this  component  is  decoupled  from  the 
axial  and  radial  components  which  are  coupled.  So,  the  tangential  direction 
is  a truly  homogeneous  direction. 

Figure  1 shows  the  Eulerian  correlation  function  for  the  tangential  com- 
ponent, denoted  by  the  solid  line,  at  r/R  = 0.4.  Least-Square-Fitting  this 
result  to  eq.  (9),  using  rc  and  tv  as  fitting  parameters,  gives  the  dashed  line 
as  shown  in  figure  1.  This  indicates  that  the  assumption  of  two  separate 
time  scales  is  justified  and  that  this  simple  model  is  capable  to  describe 
the  Eulerian  statistics.  The  correlation  function  from  the  DNS,  see  figure 
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1,  is  Fourier  transformed  to  obtain  the  velocity  spectrum,  Sv( u>).  For  the 
model  eq.  (8)  can  be  used.  Figure  1 shows  the  results  in  log-log  scale  and  in 
figure  2 a linear  scale  is  used.  Multiplying  these  spectra  with  w2  gives  the 
acceleration  spectrum  as  indicated  in  figure  2.  The  reciprocal  calculated 
values  for  rv  and  rc  are  also  shown  (arrows).  The  time  scales  are  separated 
by  approximately  a factor  10.  In  the  limiting  process,  where  tv  — > 0,  the 
spectrum  for  the  acceleration  goes  to  a constant  value  (dashdot).  This  is 
the  true  white-noise  approximation.  The  actual  DNS  results  are  far  from 
this  limiting  process.  Still,  by  using  a relatively  simply  model,  a quantifica- 
tion of  the  2 time-scales  is  possible.  In  the  near  future  we  will  analyze  the 
remaining  radial  and  axial  directions  in  the  same  way,  in  order  to  obtain 
all  (damping)  coefficients  for  the  3D  Langevin  equation.  The  result  will  be 
an  accurate  particle  dispersion  model  for  inhomogeneous  turbulent  flow.  It 
also  enables  us  to  calculate  the  universal  Kolmogorov  constant  Co- 

5.  Conclusions 

A spectral  DNS  code  is  developed  to  calculate  the  turbulent  flow  in  a 
cylindrical  geometry.  Calculated  averaged  quantities  are  in  good  agreement 
with  DNS  and  experimental  results  by  others.  Eulerian  space-time  velocity 
correlations,  calculated  with  this  code,  can  be  modeled  with  a Langevin 
equation.  In  this  way  it  is  possible  to  derive  the  damping  coefficients  for 
this  equation.  Due  to  the  relatively  moderate  Re-number  employed,  the 
white-noise  approximation  is  not  obtained. 
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Figure  1.  Left:  The  Eulerian  Space-Time  correlation  function  for  the  tangential  velocity 
component;  solid:  DNS,  dashed:  Model.  Right:  Velocity  spectrum;  solid:  DNS,  dashed: 
Model.  In  both  figures  the  arrows  represent  the  calculated  time-scales. 


Figure  2.  Left:  Velocity  spectrum;  solid:  DNS,  dashed:  Model.  Right:  Acceleration  spec- 
trum; solid:  DNS,  dashed:  Model,  dashdot:  white-noise  approximation.  In  both  figures 
the  arrows  represent  the  calculated  time-scales. 


